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ABSTRACT 

The properties of the probability distribution function of the cosmological continuous density field are 
studied. We present further developments and compare dynamically motivated methods to derive the PDF. 
One of them is based on the Zel'dovich approximation (ZA). We extend this method for arbitrary initial 
conditions, regardless of whether they are Gaussian or not. The other approach is based on perturbation 
theory with Gaussian initial fluctuations. We include the smoothing effects in the PDFs. 

We examine the relationships between the shapes of the PDFs and the moments. It is found that formally 
there are no moments in the ZA, but a way to resolve this issue is proposed, based on the regularization of 
integrals. A closed form for the generating function of the moments in the ZA is also presented, including 
the smoothing effects. We suggest the methods to build PDFs out of the whole series of the moments, or 
out of a limited number of moments - the Edgeworth expansion. The last approach gives us an alternative 
method to evaluate the skewness and kurtosis by measuring the PDF around its peak. We note a general 
connection between the generating function of moments for small r.m.s a and the non-linear evolution of the 
overdensc spherical fluctuation in the dynamical models. 

All these approaches have been applied in ID case where the ZA is exact, and simple analytical results 
are obtained. It allows us to study in details how these methods are related to each other. The 3D case is 
analyzed in the same manner and we found a mutual agreement in the PDFs derived by different methods 
in the the quasi-linear regime. Numerical CDM simulation was used to validate the accuracy of considered 
approximations. We explain the successful log-normal fit of the PDF from that simulation at moderate a as 
mere fortune, but not as a universal form of density PDF in general. 

Subject headings: cosmology: theory — dark matter — galaxies: clustering 
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1. INTRODUCTION 

One of the goals of the study of the nonlinear gravitational dynamics in an expanding universe is the 
determination of the statistical properties of the various cosmic fields that can be used to describe the matter 
distribution and motion. In the linear regime, for Gaussian initial conditions, this description is quite easy 
since it reduces to the behavior of the two-point correlation function, or cquivalcntly to the shape of the power 
spectrum. Once nonlinear effects arc taken into account this is no longer the case and the mathematical 
description of the statistical properties is more complicated. One can basically distinguish two approaches 
to describe the statistics of the non-linearities. 

The description of statistics in terms of the correlation functions has been investigated since '70s. In 
principle, a full knowledge of the matter field can be obtained from the shape of the p-point correlation 
functions (see Peebles 1980 for references). These functions are the solution of a hierarchy of dynamical 
equations, the BBGKY hierarchy. The spatial correlation functions have been measured in galaxy catalogues, 
or in numerical simulations. Progress was made to find a few low order correlation functions, theoretically, 
numerically and observationally (see Peebles 1993 for references). However, dynamical BBGKY equations 
have never been solved in general. 

Another approach which drew much attention recently is based on the probability distribution function 
(PDF) of a cosmic field at a given point or PDFs at various points, or joint PDFs of several cosmic fields. The 
one-point density PDF, P(p) dp, is the primary object of study in this approach. The discrete analogy of the 
one point PDF is given by counts in cell probabilities that is basically what is obtained after smoothing of 
the discrete field by a sharp filter. It is important to note that the PDFs contain more statistical information 
than a few lower order correlation functions. Actually the moments of the PDF are the spatial averages 
(with the same window functions) of the correlation functions. 

In practice, galaxy PDFs have been measured in various catalogues by Hamilton (1985) and Alimi, 
Blanchard & Schaeffer (1990) in the CfA survey, Bouchet et al. (1993) Kofman et al. (1994) in the IRAS 
surveys, Maurogordato, Schaeffer & Da Costa (1992) in the SSRS survey, Gaztanaga & Yokoyama (1993) 
both in the SSRS and CfA survey. The density PDF manifests significant non-Gaussian features in non-linear 
and even in quasi-linear regimes. One central question theories have to address is to determine what part 
of this non-Gaussianity came from non- linear dynamics, what part came from possible non-Gaussian initial 
conditions, and what part came from the galaxy biasing. 

To see what kind of density PDFs emerge from gravitational dynamics, investigations have been made 
in numerical simulations by Bouchet, Schaeffer & Davis (1991), Weinberg and Cole (1992), Bouchet & 
Hernquist (1993), Juszkicwicz et al. (1993), Kofman et al. (1994). 

In this paper we concentrate on the theoretical basis for the derivation of the cosmic density PDF. There 
were phenomenological attempts in the literature to design P(p) - see, for instance, Saslaw and Hamilton 
(1984); Coles and Jones (1991). Here we discuss the derivation of the density PDF from gravitational 
dynamics, and its comparison against numerical simulation. It is a difficult problem to derive P{p) for the 
general case. However, it can be studied in different regimes and approximations. One can distinguish 
different successive stages of the non-linear gravitational dynamics: quasi-linear regime, non-linear regime, 
and highly non-linear regime. The quasi-linear (or mildly non-linear) regime when a < 1 can be investigated 
by the mean of the perturbation theory, unlike other regimes. In the non-linear regime a > 1 the complexity 
of the dynamics makes analytical studies virtually impossible. Highly non-linear regime when a 3> 1 might 
take place when the hypothetic hierarchical ansatz for the correlation functions is working. Indeed, another 
feature of the gravitational dynamics - self-similar solution (Davis & Peebles 1977)- is expected to be reached 
in this regime. The statistics of this regime has been a subject of consideration in many papers. In particular 
Balian & Schaeffer (1989) give relations of crucial interest that relate the properties of one-point density 
PDF to the ones of the correlation functions. 

In this paper, we rather concentrate on the statistics in the quasi-linear regime. Fortunately, a fair 
fraction of the observational data does correspond to such a regime when the galaxy surveys are smoothed 
with a large enough radius (say £ 8/i _1 Mpc). In the quasi-linear regime the two-point correlation function 
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is mainly determined by the initial conditions. The manifestation of the nonlinear features will be seen in the 
higher order correlation functions, or equivalently in the departure of the shape of the distribution functions 
from Gaussian distributions. We have now various tools at our disposal to study the mildly nonlinear regime 
in details, so that specific predictions can be made that are directly derived from first principles. 

One of the most successful approximations to describe the early nonlinear dynamics is the Zel'dovich 
approximation (Zel'dovich 1970). In this approximation, the displacement of the particles is extrapolated 
from its behavior in the linear regime. The whole local statistical properties of the field at the position of a 
particle is then determined by the statistical properties of the initial displacement field. For Gaussian initial 
condition, it is possible, for instance, to compute the density PDF by a change of variable starting with the 
initial statistical properties of the displacement. Actually, as it will be recalled in this paper, a lot more 
statistical properties can be derived with such an approximation (Kofman 1991, Kofman et al. 1994). 

The second method is based on the calculations of the large-scale cumulants of the cosmic field from 
perturbation theory. Indeed, when one assumes Gaussian initial conditions, it turns out that it is possible to 
derive the leading term (with respect to the rms density fluctuation) of each cumulant (Bcrnardeau 1992). 
These first order contributions are expected to dominate the exact value of the cumulants as long as the rms 
fluctuation is accurately given by the linear theory. The shape of the density PDF is then obtained by a 
reconstruction of this function from the generating function of the moments. 

The mathematical forms of the density PDFs derived analytically in these two methods are quite 
different. We compare the forms of the PDFs derived in the quasi-linear regime earlier in our papers, 
and show its mutual agreement. As this subject is in rapid development, many questions and controversies 
on the density PDF and moments are accumulating in the literature. Therefore we try to clarify many points 
related to the PDFs, generating functions, moments and smoothing in the quasi-linear regime. Presumably, 
there is no simple universal formula for P{p) in general. However, one of the practical outputs of our paper 
is to provide a justification for the use a log-normal distribution as a successful fit for the actual cosmic 
density PDF. This fit is only accurate for the cosmological models based on the Gaussian initial statistics 
and realistic power spectrum, but not in the general case. Another practical output is a new view on the 
measurement of the lowest moments - skewness and kurtosis. Common belief is that they are affected by the 
high density tail which is difficult to measure. We will demonstrate that these moments can be evaluated 
from the form of the PDF around its maximum. 

We will use the concepts of moments, cumulants, PDFs and generating function throughout the paper. 
The formal definitions of all of them are collected in Appendix A. However, we must stress that the 
introduction of quantities such as the generating functions of moments and cumulants is not an extra 
mathematical exercise but an useful tool to find the measurable statistics of the cosmic fields. There is a deep 
connection between generating function and the non-linear evolution of the spherical overdense fluctuation 
in the dynamical model of the gravitational instability. Hence, in most cases the generating function can be 
derived directly from the dynamical equations. We will present the closed forms of the generating functions 
in different approximations in the quasi-linear regime. 

Throughout this analysis we got a series of new results: the generating functions of cumulants and 
finite moments in the Zel'dovich approximation, both with and without final smoothing; two methods of 
reconstruction of cosmological density PDF from, correspondingly, the whole (Laplace transformation) and 
partial (Edgeworth decomposition) series of moments; systematic mutual comparison of all these methods 
and their comparison against the calculations from N-body simulation; the range of validity of their fitting by 
the log-normal distribution; the extension of the derivation of density PDF in the Zel'dovich approximation 
for an arbitrary initial condition. 

Section 2 of the paper is devoted to the ID dynamics, in which the mathematical content is simpler 
than for the 3D dynamics. We apply general methods and mathematical tools in this case. Section 3 is 
devoted to the statistical properties that can be derived from the use of the Zel'dovich approximation for the 
3D dynamics. In Section 4, we give the properties of the exact dynamics derived from direct perturbative 
calculation in the single stream regime. Section 5 is devoted to comparisons with numerical simulations. In 
the conclusive Section 6 we summarize the approximations that have been made and the results in a detailed 
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2. AN ILLUSTRATIVE EXAMPLE: THE ID DYNAMICS 

In this part we aim to present the analytical tools developed in the mildly nonlinear regime in case of 
a simpler dynamics. We assume in this part that the fluctuations are only one-dimensional. The general 
methods suggested here will be used in the next sections in the 3D case. 

2.1. The construction of the density PDF from the Zel'dovich approximation 

It is convenient to use the equation of motion in the Lagrangian description. In such a case the dynamics 
is described by the displacement, S(q, t), of each particle from its initial position q. Its current Eulerian 
comoving position, x, is then given by 

x = q+S(q,t). (1) 

For ID case the mildly non-linear dynamics is quite simple. The reason is that the force exerted by a density 
perturbation over a given particle is independent of its distance to the particle. Therefore, before any shell 
crossing, the displacement of each particle depends on its Lagrangian position q only. Then the displacement 
field for the growing mode can be factorized 

S(q,t) = *(q)D(t). (2) 

The 3D generalization of this form of the displacement is the Zel'dovich approximation (ZA). In ID case the 
ZA is then identical to the exact dynamics. In relation (2) the displacement factorizes in a spatial function, 
*(q), which depends on the initial conditions, and a universal time dependent function D(t) which contains 
by definition the time dependence of the growing modes. Let us call A the derivative of — W(q) with respect 
to q. The local comoving density (in units of the mean density, g = p/p) is then, by the mass conservation 
equation, given by 

M - idscr (3) 

The reconstruction of the density PDF then is based on the statistical properties of A . For Gaussian 
initial conditions, Ao simply obeys Gaussian statistics. The derivation of the shape of the density PDF is 
then just a matter of change of variable, and as it can be seen in (3) the density contrast will not be normally 
distributed. This is general feature due to the nature of this quantity. The positive density contrast can reach 
any large value while the negative density contrast is restricted by g > 0. Hence the probability function 
P(g, D(t)), meaning the fraction of volume with a given value of density, is expected to be very non-Gaussian 
even in the quasilinear stage. However, note that a simple change of variable leads to the density PDF at a 
given point, in Lagrangian space. The distribution in Eulerian space should take into account the fact that 
a given amount of matter in a dense spot occupies a small volume. The density PDF in Eulerian space is 
then obtained by divided the density distribution in Lagrangian space by the density, and multiplied by the 
numerical factor which controls the normalization. This factor is related to the number of streams (N s ), in 
cases (7^1 under consideration it is very close to unity, and we ignore it (e.g. Gurbatov et al. 1991). 



Then the density PDF in ID reads 



Pi D (g)dg = 



1 



(2tt ( t 2 ) 1 /2 



exp(-A_)+exp(-^_) 



A = £>A = 1--, X' = D\'=1 + -. 

g g 



dg 

(4) 



Note that o = gq D(t), where do is the variance of the initial linear density contrast. In the limit of small o 
the distribution (4) transits to the Gaussian form. The high density asymptota g~ 3 that is seen in equation 
(4) is induced by the caustics. 
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2.2. Density PDF and moments 

As recalled in the introduction, the knowledge of the shape of the PDF and the moments are intimately 
related. The very example we are presenting now is here to point out that this is not completely true. Indeed, 
none of the moments of the previous distribution are finite, because of the £) _3 -asymptota! It can be easily 
checked that even for an arbitrary small a the density PDF behaves like -\/2/7r/crexp(— \/2a 2 )g^ z at high 
density. However, physical density in caustics is finite. We expect that a physical process will operate as 
sort of regularization of the high density tail of P(g). For example, Zel'dovich & Shandarin (1982) showed 
how physical properties of the hot dark matter lead to the finite density in the caustics. Obviously we 
do not know what will be the exact form of this regularization in general case, so we simply describe the 
regularization process by a cut-off at large density (where the physical regularizing processes are thought 
to occur) that makes the moments finite. The shape of this cut-off is somewhat arbitrary, but we will see 
that, to some extent, a fair fraction of the properties of the moments do not depend on the procedure that 
has been adopted. We chose two toy models of distribution, with the sharp cutoff and with the exponential 
cutoff, in the following form: 

Piegi(Q) dg = -^-Pid(q) dg if g < g c , 
^rcgi(e) dg = if q > g c , 

and 

-Pro g 2(£>) dg 

X 

The coefficients C\ and C2 arc simply normalization factors. The adopted values of g c will be quite large 
(above 5) insuring that the shape of the density PDF is not changed in the domain of interest that is for 
< g < several. The coefficients C\ and C2 are very close to unity in these cases. In Fig. 1 we present the 
distribution (4) and regularized distributions (6). 

Once such a regularization is made, the moments of the distribution functions (S p ) = J °° dg P(g)(g— l) p 
can be easily calculated numerically. We present the results in terms of S p parameters adopted in the 
literature. These parameters are defined through the cumulants (see Appendix A, Eq. [A2]): 

S P = (8 p ) c /o 2 <*- 1 \ (7) 

This is a generic definition which is not based on any a priori assumption. The S p parameters would be 
constants in the hierarchical ansatz, but in general they can be seen as functions of a. 

The parameter S3 multiplied by a is the skewness, and S4 multiplied by a 2 is the kurtosis. The numerical 
calculations of the parameters 53 and 5*4 for the regularized distributions (5) and (6) for the exact Zel'dovich 
solution is shown in Fig. 2. The measurements have been made by calculating the actual second, third, 
and fourth moment of the regularized density PDF, and by taking the appropriate ratios. We see that these 
coefficients are finite. The shapes of 53(a) and 84(a) are universal in the interval of a from up to the value 
~ 0.2. These shapes are independent of the form of the cutoff, as well as on the parameter of truncation. In 
the limit of a — > we found S3 w 6, and S4 w 72. In the next section we will confirm these figures and their 
universality by analytical calculations. We conclude that the found values of 83(a) and 84(a) for a ^ 0.2 
are the proper moments in the ID case. There is also some universality of these coefficients for a ^ 1.2, in 
the sense that it does not depend on the parameter g c for each panel. However, it does depend on the form 
of regularization. We suggest that the found values of S3 (a) and 84(a) for a £ 1.2 do not reflect the proper 
moments, but rather the shape of the adopted cutoff and thus do not have any physical meaning. 

2.3. Calculation of S p for small a 

It is possible to calculate the S p series analytically when the variance a of the distribution function is 
small. This calculation requires the introduction of the generating functions of moments and cumulants (see 



(5) 



= i- PMg) 
= g- g c - 



1 + 



cxp(a;) 



cxp(a;) + exp(— x) 



(exp(-g/g c ) - 1) 



dg, 



(6) 
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Appcndix A). The formal definition of the generating function of the moments (5 P ) is 

M(h) = 1 + J2{5 p )^t (8) 
P =i p - 

and the generating function of the cumulants (S p ) c is 

Here ^ is an auxiliary parameter. Using its definition, one can relate M.{li) to the shape of the density PDF: 

/>oo 

M(n) = expC(/i) = / deP(e)exp([ e -l]/i). (10) 
Jo 

Since we know the density PDF in the quasi-linear regime in the ID problem, in principle, we can obtain 
the generating function C(/x) substituting the formula (4) into integral (10). When a <C 1 we can integrate 
(10) in a closed form, using the steepest descent method. The saddle point of the exponent is given by the 
equation 

1 d (\{gf 



a 2 dg 



+ fi = (11) 



with X(g) = 1 — l/g. Note that the solution of this equation, A s , is such that DX S is the function of the 
combination /xcr 2 only. In ID case the solution of equation (11) can actually be obtained straightforwardly. 
We will treat, however, this equation in a different manner, which will be much more appropriate in 3D case. 
We introduce the function Q(r), 

g(r) = -^ 1, so that g = Q{-\) + l. (12) 

1 + T 



After substitution of formula (12) into equation (11), the solution of equation (11), A s , is then given by 
the implicit equation 



2 d<? 



D\ s = — 
At 



(13) 

t = -D\ s 



The generating function C(/x) of the cumulants is obtained by taking the logarithm of the integral (10). In 
the low a limit it leads to retain only the term that is under the exponential at the saddle point position. 
We then obtain from (10), 

C<j*) = »g(-r s ) - ^r, t s = —D\ s (14) 

where the saddle point A s (or equivalently t s ) is given by the equation (13). 

The reason we introduced Q(t) is the following. Equations (13) and (14) have a structure of the Legendre 
transformation from Q(t) to C(/i) with the variable of the transformation equal to unity afterwards. It turns 
out that the function Q{t) is the central object derived from the basic dynamical equations in the perturbation 
theory calculations (Bernardeau 1992). In the perturbation theory Q(t) is defined as the generating function 
of the other averages - vertices (see Sec. 4). It is remarkable that it corresponds to the dynamics of the 
"spherical" collapse - one dimensional in this case: Q(— A) gives the density contrast of a collapsing object 
of linear overdensity Ao, A = D(t) Ao, as function of time. We will discuss this derivation in general case 
in Sec. 4. The reason the Legendre transformation emerges is that two generating functions Q(r) and 
C{n) considered as the sum over the corresponding averages are connected through that transformation. 
Hence, note that in the ID case expression (14) for C(^i) can be obtained from the direct perturbation scries 
expansion of the basic equations as well. 
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The crucial observation from the equation (14) is that the combination C(/x)//x is a function of the 
combination [ia 2 only. What does it imply? The term in uP , the coefficient of which is the p th cumulant we 
are looking for, is then proportional to a 2 ^ 1 ^. It rigorously demonstrates a scaling relation between the 
cumulants in the small a limit (i.e. at large scales). The expression of these coefficients can be obtained 
with an expansion of the function C(/x) with /x: C{fi) = n. 2 cr 2 /2 + fj, 3 a 4 + 3u, 4 <j e + ... . Then, for the skewness 
and kurtosis, for instance, we get 

S 3 (0) = 6, (15a) 

and 

5 4 (0) = 72. (156) 

These results can be extended by the further perturbative calculation that takes into account the next 
terms of the expansion in the steepest descent method. It is then possible to get the first a corrections for 
the expression of the skewness and the kurtosis. We then get 

S 3 (a) =6 + 24o- 2 + 0(cr 4 ), (16a) 

and 

S 4 (a) = 72 + 810cr 2 + O(cr 4 ). (166) 

These expansions are not affected by the shape and parameters of the cut-offs. Theoretical curves (16) 
derived for small a are plotted on Fig. 2, and indeed osculate the universal numerical curves for small a. 

One of the challenges of the study of the large-scale structure formation is to find the a dependence 
of the S p parameters. From the theoretical point of view, it is not clear whether the a dependence can 
be accurately determined with perturbation theory in the single stream approximation. Multistreaming 
may change the behavior in an unknown way. But we present here the results in ID case to stress that 
perturbation theory does not prove at all that the parameters S p are constant when a < 1. The important 
property is that, for Gaussian initial fluctuations, they have finite limits at small a. Note that the expansion 
of S p as a function of a similar to formula (16) has never been done for the 3D case. 

2.4. The reconstruction method from the generating function of the cumulants 

In the previous section we solved the problem how to find the moments from the known PDF. When 
the whole series of the cumulants is known it is possible to to construct the density PDF from the generating 
function of the cumulants by inversion of the relation (10) (the inverse Laplace transformation): 

j / i +ioo 

P(e)dQ = ^ / d M exp {-C(p) - {q - 1) M ] , (17) 

27T1 J-\oo 

where the integration is made in the complex plane. 

In practice we can find the moments and generating function in the limit of small a. For the further 
analysis let us define the function ip(y) by the relation 



00 (-\\p- 1 



<p(y,a) = > S p (a)^—yP = y - a 2 C(-y/a% (18) 



Pi 

P =i 



where we define Si = 1, S2 = 1. The advantage of this new function is that it is finite for an arbitrary y in 
the limit of small a. This function is the generating function of the parameters S p (a). Then from (17) and 
(18) the density PDF is given by 



2-irr-' 



J —loo 



(19) 
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Balian & Schacffer (1989) used this form for the hierarchical ansatz, when there is no dependence on a in 
tp(y,cr). In the small a limit the a dependence that may be contained in tp(y,a) vanishes, ip(y,cr) — > <p(y), 
and the results of Balian and Schaeffer (1989) can then be used here, too. 

Note that general formula (19) does not assume, a priori, that a is small (as for the Edgeworth expansion 
in the next subsection). It supposes, however, that the function ip(y,a) that is used is valid in the domain 
of application. The reconstruction formula is of obvious interest when the function ip(y, a) can be derived 
from the first principles as it has been shown when a = in ID case. We will see that it is also the case for 
small a limit in 3D case. The reconstruction formula (19) is thus of general interest. 

In Fig. 3 the density PDF in ID case is calculated numerically with the reconstruction formula (19), 
where ip(y) is calculated from (14) and (18) in the small a approximation, i.e. ignoring cr-dependence 
in ip(y,<r). It is then compared to the original shape of the density PDF (Eq. [4]). We find that the 
reconstruction method based on the generating function works very well for the density interval < g 2, 
slightly worsening as <r increases. The shape of the high-density tail at g ^ 2 of the PDF from the 
reconstruction method differs from that of the actual PDF. For the form (19), the density PDF has an 
exponential cutoff as exp[— £>/(x*cr 2 )] with x* = 27/4 = 6.75 (see Bernardeau 1992), which is very different 
from the g~ 3 cutoff of formula (4). This discrepancy is just due to the fact that we ignore the cr-dependence 
in the parameters S p . Therefore the discrepancy is increasing with a. It is quite interesting that for moderate 
a, the cr-dependence affects the high-density tail of PDF, but does not affect its shape around the maximum, 
nor the low-density tail (see the Table in Sec. 6). 



2.5. Reconstruction of PDF from a few cumulants: The Edgeworth expansion 

In the previous section it has been shown how it is possible to derive the shape of the density PDF out of 
the generating function. In this section we report a method that can be used to recover the shape of the PDF 
when only a limited number of cumulants is known. In practice we may have only a few lowest cumulants, 
such as the skewness and the kurtosis. In the case of the weakly non-linear dynamics, when slight departure 
from the initial Gaussian distribution is expected, one can use the general decomposition series around the 
Gaussian PDF induced by the first non-zero cumulants. This decomposition is known as the Gram-Charlier 
series (Kendall & Stuart, 1958). Longuct-Higgins (1963) applied the Edgeworth form of the Gram-Charlier 
series to the statistics of the weakly non-linearities in the theory of 2D sea waves. Inspired by this paper, we 
suggest to use the Edgeworth's decomposition for the ID and 3D cosmological density PDF (as reported in 
Kofman 1993). We understand that similar ideas were independently suggested by Juszkiewicz et al. (1993). 

The Edgeworth expansion can be derived from the form (19) of the density PDF. Assuming that the 
density contrast 5 is of the order of cr and small, the relevant values of y are also of the order of cr and are 
thus assumed to be small. It is then legitimate to expand the function ip(y) 

v{y, a) « y \v 2 + § y 3 § y 4 + § y 5 ± . . • , (20) 

where S p = S p (<r). To calculate the density PDF, we substitute the expansion (20) into the the integral (19). 
Then we make a further expansion of the non-Gaussian part of the factor exp [— ip(y)/a 2 ] with respect to 
both y and cr assuming they are of the same order, see Appendix B for details. 

Finally we obtain the so-called Edgeworth form of the Gram-Charlier series for density PDF 



1 + a^H 3 (u) + cr 2 ( §ff 4 („) + $H 6 („) 



1 120 W 144 K ' 1296 V ' 



6 w V 24 w 72 
dS, 



(21) 



where H n (v) are the Hermite polynoms (see Appendix B), v = S/a. This is a universal form for any 
slightly non-Gaussian dynamical models, i.e. when a is small and S p are finite. The actual forms of the 
parameters S p = S p (a) which depend on particular dynamics, affect the expansion (21) with respect to a. 
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Longuct-Higgins (1963), Juszkiewicz ct al. (1993) and Kofman (1993) gave it up to cr 2 -correction. We give 
the Edgeworth expansion up to cr 3 -correction, for which the CT-dependence of S3 (see, for instance, eq.[16a]) 
has to be taken into account. The resulting approximate form (21) is a Gaussian distribution multiplied by 
a corrective function - in the square brackets - expanded with respect to ct. The zero order of this expansion 
gives simply a Gaussian distribution; the first order corrects it by taking into account the skewness; the 
second order by taking into account both the skewness and the kurtosis, etc. 

Thus, it is possible to get an approximate form of the density PDF from a few lowest known cumulants. 
This method is irreplaceable when only a few cumulants have been derived from the first principles. However, 
it is important to note that this expansion has been made possible only for slightly non-Gaussian regime. 
The validity domain of the form (21) is limited to finite values of 6/ a. 

Let us now apply the Edgeworth expansion (21) to the ID gravitational dynamics. We have to take 
into account the ct dependence of 5*3 given in equation (16a), and then obtain, 



We plot the density PDF in ID case reconstructed from the Edgeworth expansion (22) in Fig. 4, and 
compare it with the actual PDF (4), for ct = 0.1, and ct = 0.3. We can see that a few iterations of the 
expansion (22) reproduce the peak of P(5) in the interval \S\ £ 0.5 around it for small ct relatively well. It 
reproduces well the shift of the maximum towards the low density. It completely fails to reproduce P(S) 
outside of this interval where spurious wiggles appear. For a given value of ct, each next ct iteration quite 
slowly improves the approximation. Unfortunately, the method is rapidly worsening as ct increases, and in 
practice is useless for ct £ 0.5. 

The usual measurements of the lowest cumulants are significantly affected by the high density tail of 
the PDF, i.e. the rare events. It is interesting to note, that in the context of the reconstruction methods, 
the lowest cumulants alone are responsible for the shift of the peak of P(S). Therefore the measurement of 
the shape of the PDF maximum, which statistically is more robust, can provide an alternative method of 
evaluation of the lowest cumulants. 



We consider here the statistics of the cosmic density field in 3D Zel'dovich approximation in a similar 
manner as it was done in ID case. The most important change is that the Zel'dovich solution no longer 
reproduces the exact dynamics but is an approximation. It is, however, thought to be a good description, so 
that it is worth investigating the statistical properties that can be obtained in this approximation. We report 
here a different derivation of what have been done by Kofman (1991) and Kofman ct al. (1994) for Gaussian 
initial fluctuations. The new method (Kofman 1994) also allows to extend the results to non-Gaussian 
adiabatic initial fluctuations. 

3.1. Zel'dovich approximation for filtered initial fluctuations 

For the sake of simplicity we assume that the universe is filled by the pressureless matter. The growing 
mode of adiabatic perturbations is D(t). Let x, v = a dx/dt, p(t, x) and <p(t, x) be, respectively, the comoving 
coordinates, peculiar velocity, density of dark matter and peculiar Newtonian gravitational potential. It is 
convenient to introduce a new time variable, D(t), then a comoving velocity is u = dx/dD. The growing 
mode is non-rotational, so that the velocity field is potential. Let $ be the velocity potential so that u = V x $. 



P(S)d5 




(22) 



3. THE 3D DYNAMICS WITH THE ZEL'DOVICH APPROXIMATION 



The gravitational dynamics of the cosmological system is complicated and requires the N-body 
simulations. For interesting cosmological models such as the CDM scenario, the structure formation looks 
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like complicated hierarchical pancaking and clustering from very small to large cosmological scales (e.g., see 
Shandarin & Zel'dovich 1989, Kofman et al. 1992). However the gravitational clustering at sufficiently large 
scales R can be considered in the quasilinear theory in a single stream regime ignoring small scale details. 
For this goal we use the Zel'dovich approximation but apply it for the filtered initial gravitational potential. 
This approach, sometimes called the truncated Zel'dovich approximation, was used for different purposes in 
papers Bond and Couchman (1986), Kofman (1991), Shandarin (1992), Kofman ct al. (1992), Coles ct al. 
(1993), Kofman et al. (1994). From a mathematical point of view the truncated Zel'dovich approximation 
means just that the Zel'dovich approximation is applied to the truncated initial potential to ensure being in 
the single stream quasilinear regime. 

To describe the motion of particles we can introduce the tensor of the velocity derivatives Sij(x,t) — 
— V Xi Vj in the Eulerian space; for the potential motion it is reduced to Sij — — V Xi V Xj $. Let Xi be 
its eigenvalues. The field of the SV, (i)-tensor evolves in time, its initial value (in the Lagrangian space) 
coincides with the Lagrangian deformation tensor Dij = — V qi V qj $o- From dynamical equations one can 
easily show that the eigenvalues Xi(t) are related to their initial values Xoi, by Xi(t) — Aoi/(l + D(t)Xoi). 
The local density can be obtained by the inverse of the Jacobian of the transformation between q and x, so 
that 

Q= \(1 - DX al )(l - DX 02 )(l - DX 03 )\- (23) 

In the ZA the statistics of the evolved field can then entirely be obtained by the statistics of the initial local 
density g and the initial eigenvalues A -s. For adiabatical perturbations go = 1. 

3.2. Joint PDF in the ZA for arbitrary initial statistics 

In the last section 3.1 all relationships were obtained without making any assumption on the initial 
statistics. The cosmic density PDF can then be obtained from the initial joint PDF of all involved cosmic 
fields: Wo(qo, Aoi, A02, A03, uo, $o) d^o dAoi dAo2 dAo3 d 3 uo d$o- That the statistics can be completely made 
from the statistical behavior of all variables is entirely due to the use of the Zel'dovich approximation. If 
this approximation were released it would no longer be true. 

The density PDF can be obtained in general case of arbitrary initial condition by integrating the 
combination £[(?|(1 — DAoi)(l — DA02XI — D\os) \ — l] x Wo(qo, Aoi, A02, A03, Uo, $0) over all involved variables 
except density. This is one of the new results of our paper. Density PDF for non-Gaussian initial fluctuations 
will be considered in a separate paper (Kofman 1994). 

3.3. Density PDF for Gaussian initial fluctuations 

In this Section we study the statistics of the continuous cosmological fields evolving from the initial 
Gaussian fluctuations, which is the general frame for most cosmological models. For the Gaussian initial 
conditions we can omit uo and $0 m the initial joint PDF and write it as 

W (go, Xqi, X02, X 03 )dg dA i dA 2 dA 03 = S(g Q - 1) dg M (A i, A 02 , A03) dA i dA 02 dA 03 . (24) 

The first factor - the Dirac <5-function - is the initial density distribution function, which corresponds to 
the perfectly homogeneous density distribution g = 1. This is just the formal limit of the Gaussian density 
distribution with a — > 0. This factorization is expected to take place for cosmological models with small 
adiabatic initial fluctuations. The second factor is the joint distribution function of the eigenvalues of the 
initial deformation tensor for an initial Gaussian displacement field (Doroshkcvich 1970) 

dAoi dAo2 dAo3 Mo(Aoi, A02, A03) = dAoi dAo2 dAo3 
5 5 / 2 27 

x — — g- (A i - A 02 )(A i - A 03 )(A 2 - A 03 ) exp 



-I 



(12 



(25) 



where J i = A i + A 2 + A03 and J 2 = A iA 2 + A iA 3 + A 2A 3- 
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The shape of the density PDF can then be obtained by the change of variable g to g and by integrating 
over Ao-s in (24). We then have 

P(q, D)dg = dg f d\ 01 dA 02 dA 03 6[g\(l - DX 01 )(1 - D\ 02 )(l - DX 03 )\ - l] 

J (26) 

x M (A i, A 02 , A 03 ). 

We get the formal expression for the Eulerian density PDF in the ZA. 

Substituting the expression (25) into the integral (26), after some tedious algebra, we can reduce the 
integral (26) to a simpler one-dimensional integral which has to be evaluated numerically, 



P(g,D)dg 



_ 9-5 3 / 2 dg 
"47r7V s(? 3 cr 4 y 3 (|)V3 



ds 



-(s-3) 2 /2a 2 



^ (g-/3?/2a 2 + e -/? 2 2 /2a 2 



-/3 3 2 /2 



(27) 



COS 



2 . N 1 

— {n — l)ir + — arccos 



54g 3 
gs 3 



1 



where the parameter a(t) = D(t)ao is the standard deviation of the density fluctuations g/g in the linear 
theory, and N s is the mean number of streams, N s = 1 in the single stream regime. The expression (27) was 
derived by a different method earlier on by Kofman (1991), Kofman et al. (1994). 

This 3D formula corresponds to the formula (4) derived in ID. The analytical expression is obviously 
different, but the method to build it follows the same scheme: the density PDF is obtained from the joint 
PDF of the eigenvalues Ao of the initial local deformation tensor. In the limit of very small a the formula 
(27) is reduced to the Gaussian distribution. We plot the PDF calculated numerically from formula (27) in 
Fig. 5. For moderate values of a density PDF calculated in the ZA is in good agreement with the PDF from 
N-body CDM simulations, sec Kofman et al. (1994), and also Fig. 7. One can expect formula (27) works 
even better for those models, for which the pancaking is more pronounced. 



3.4. Calculation of S p for small a in the ZA 

The 3D density PDF in the ZA also has the caustic-induced g _3 -asymptota. As a result, the moments 
cannot be formally defined in the ZA for any given value of a. However, we can regularize the PDF in the 
Zel'dovich approximation by cutting off the high density tail as we did in Sec. 2.2. for ID case. Then we are 
able to define the cumulants and derive their large scale properties. We use the same two regularizations, 
the sharp cutoff and the exponential cutoff as for the ID case, see equation (4). The parameters S3 and S4, 
defined as in relation (7) and calculated with numerical integrations of the moments, are plotted in Fig. 6 
as functions of a. They exhibit a very similar qualitative behavior than for the ID case. Numerical shapes 
of S p (a) are universal for small a 0.3, i.e. independent of the form of the cutoff and on the parameter of 
truncation. However, their quantities in the small a limit are changed compared with ID case. In 3D case 
we found S3 « 4 and S4 w 30. For er > 0.3, S3 and S4 depend on the shape of the regularization functions, 
therefore the moments in the ZA for moderate and large a are poor-defined. 

In the small a limit more advanced analytical progress can be done for the Gaussian initial fluctuation. 
Grinstcin & Wise (1987), Munshi & Starobinsky (1993) calculated S 3 (0) = 4; Bernardeau et al. (1994), 
Catelan & Moscardini (1993) calculated 5*3(0) and ^(O) = 272/9, all from the perturbation theory around 
the ZA. The new result we present in the rest of this section is the analytical derivation of the generating 
function of the cumulants, and consequently, in principle, all the cumulants themselves, in the ZA. 

The most straightforward method of the calculation of the generating function C (p) and cumulants is 
the same as for the ID case. Let us use the integral (10) for C(n), where for the density PDF in the ZA we 
use the integral (26), which includes the joint PDF of the eigenvalues Aoi and the S- function of the density. 



- 13 - 



Integrating the ^-function over the density, we get the integral over Aoi-s 



5 5 / 2 27 r dA 01 dA 02 dA 03 

87ra§ J |(1 - £>A i)(l - £>Ao2)(l - D\ 03 )\ 



eXpC(fl) = c _ g / — — — — — — — — (A i - A 02 )(A i - A 03 )(A 2 - A 03 ) 



x exp 



(28) 

J_(„ T 2 _ 15 A ( 1 

- ^ ^3J 01 2 J02 J + M ^ |(1 _ DXoi){1 _ DAo2)(1 _ DXo3)l 



In the limit of small a we can apply the steepest descent method to evaluate this integral. This is a triple 
integral that, however, can be calculated since Aoi-s are involved in the integrand in symmetric combinations 
only. There are three equations for the saddle point (A s i, A s2 , A S 3) in the 3D A^-space, which admit a simple 
symmetric solution 

A01 = A02 = A03 = A s , 



3 -DA, = 



a 2 M (29) 



(1-^A*) 4 ' 

Note that D\ s is a function of the combination a 2 11 only, as it was in ID case. This fifth-order equation 
cannot be solved analytically. Despite that, the introduction of the machinery of CJ(r)-function described in 
Sec. 2.3 helps to overcome the problem. 

The saddle point equation (29) can be rewritten in terms of £ z (r)-function. Indeed, the equation 

T s = ^ 2 ^t(t s ), T S EE -3L»A S , (30) 
is reduced to the algebraic equation (29), if we choose 

6ZM =arw- L <31) 

Then the leading factor of the integral (28) for small a, emerging from the steepest descent method, gives 
us the expression for C(/i): 

C z (ii)=^ z (T s )-^, (32a) 

or equivalently 

<P z (y) = y + yG z (r s ) + ?f, Ts = -y^g z (r). (326) 

The formulae (31), (32) are the 3D counterparts of (12) and (14) for the ID case. Again, equations 
(32) have the structure of the Legendre transformation from Q z ir) to C z (p) with the variable of the 
transformation equal to unity afterwards. It is remarkable, that the formulae (31) and (32) can be derived 
independently from the perturbation series based on the dynamical equations of ZA, similar to the method of 
Bernardeau (1992) based on the perturbation series of cosmological equations (see also Sec. 4), and what we 
discussed in Sec. 2. With this approach we obtain form (32) for C z (/j,) where the function Q z {t) describes 
the collapse of the density contrast of a symmetric perturbation of linear overdensity — r in the equations 
corresponding to the three dimensional ZA. The solution of this problem is given by formula (31). Note, 
that the form (31) has been derived here by a totally different method! 

Using the expression (32), we can derive analytically the full series of the S p (0) parameters in the ZA. 
For instance we have 

S z (0) = 4, (33a) 
979 

Sf(0) = -f « 30, (336) 
S §(0) = ^ ~ 342. (33c) 
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The first two figures are represented by circles in Fig. 6, and give the small a limit of the numerical curves 
S3 (a) and 84(a). The next terms of the asymptotic series of the integral (28) would allow, in principle, to 
derive analytically the a expansion of the parameters S p (a), similar to ID decompositions (16). We leave 
this exercise out of the scope of the paper. It is interesting to note that the skewness, kurtosis, etc, in the 3D 
case are systematically smaller then corresponding numbers of the ID case, c.f. (33) and (15). As numerical 
curves of Figs. 2 and 6 show, the a dependence (where it is well-defined) is also weaker in 3D case than in 
ID case. In the ZA the departure from the Gaussian distribution is than faster in ID case then in 3D case. 

3.5 Effects of the final smoothing in the ZA for small a 

In the previous section we found the parameters S v and the generating function in the limit of small a in 
the framework of the ZA. Using the machinery of the reconstruction, described in Sec. 2., we can obtain the 
density PDF in this limit from the generating function (19), or from the Edgeworth asymptotic expansion 
(21). However, it has little sense since we know the PDF (27) for an arbitrary a in the ZA. What is of greater 
interest is to use the reconstruction technics when the final smoothing effects in the generating function and 
S p are taken into account. It can be done for small a only, and allows, in principle, to reconstruct the density 
PDF in the ZA with final smoothing. 

To get the PDF for the continuous field for different a, we need to filter the Eulerian density field, 
cither in observational surveys, N-body simulations, or analytical approximations. The truncated Zel'dovich 
approximation, introduced in Sec. 3.1, is based on the smoothing of the initial fluctuations. We recall, that 
the initial smoothing is inevitable anyway, otherwise we get the multiple streaming regime, for which the ZA 
is not applicable at all. The fact that the order of the smoothing and the dynamical evolution not commute 
was noticed by Kofman et al. (1994). It is quite complicated to incorporate analytical approaches given 
in that paper or in Sec. 3, with the final smoothing. The reason is that the smoothed density does not 
only depends on the local eigenvalues of the deformation tensor at one point, but on the behavior of the 
deformation tensor in a whole smoothing large area. 

When the filtering is taken into account, the PDF is expected to slightly depend on the shape of the 
power spectrum, since it is known after Goroff et al. (1986) that the parameters S p depend on it. On the 
other hand, the PDF (27) does not have any dependence with the power spectrum and is only defined by 
the value of the rms density fluctuation a, showing some limitations for the practical use of this result. For 
the sake of simplicity we will assume in the following, that the dependence on the power spectrum can be 
reduced to the dependence on the effective index n at the scale at which the field is filtered. The index n is 
the logarithmic derivative of a with the scale R. 

There was an attempt by Padmanabhan and Subramanian (1993) to take into account the additional, 
final smoothing. They made an additional approximation within the ZA, assuming that the typical collapse 
is spherical. Unfortunately this assumption is valid for a small fraction of the Lagrangian volume, and 
the resulting distribution poorly fit the numerical PDF. Juszkiewicz et al. (1993) suggested an interesting 
phenomenological conjecture in order to take into account the final smoothing effects. They pointed out 
that the a and n dependence of the density PDF (for small a) is mainly contained in the £3 a combination. 
In the context of the ZA this property leads to the substitution a — > a 5 3 (n)/4 in the density PDF (27). 

In this Section we present the generalization (in the limit of small a) of the generating function (32) in 
the ZA that takes the final smoothing into account. The top- hat filtering is the simplest model of smoothing 
for which a complete analytical study can be done. We report the results of a general method developed by 
Bernardeau (1994) and applied here for the ZA. The basic idea is that the non-linear contributions to the 
S p parameters are not affected by the filtering at a given mass Mo scale (top-hat filtering in the Lagrangian 
space). This property is valid for the ZA as well as for the exact single stream cosmological dynamics, 
but may not be true for other dynamical approximations. The filtering effects at a given radius Rq can 
then be calculated from a transformation of the density PDF from Lagrangian space to Eulerian space: 

qPe{q 1 Rv)&Q = Pl{q, M )dg 7 where M = AttRq/3 (see Bernardeau 1994 for details). It is striking 
to see that the result in terms of a (/-function can be derived from that through (32) and the reconstruction 
formula (19). For the ZA with final smoothing this function will be denoted as Q zs . The mathematical 
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transformation to obtain Q zs from Q z (corresponding to the case without smoothing) is then given by the 
formula, 

where Rq is the filtering scale and a(Ro) is the rms density fluctuation at this scale. The generating function 
of the cumulants is then given by 

v zs (y) = y + yG zs (r) + \r\ r = -y^G zs (r) 

For instance the coefficients S zs (0) and S zs (0), for a power law spectrum of index n, are 

S zs (0)=4-(n + 3), 

The smoothing correction given by formulae (34), (35) are quite complicated, but final formulae (36) are 
rather readable. 

Fortunately, for one of the most intersting cases n = — 1, corresponding to the CDM power spectrum 
on the galaxy clustering scale, there is a simple solution of equation (34) 

G zs (t) = (l-r/3) 3 -l, (37) 

note difference with eq. (31) for Q z (t). Using this simple form, it is possible then in case n = — 1 to build 
the PDF from the reconstruction formula (19). We present the result for PDF for this interesting particular 
case in Fig. 5 and compare it to the shape of the density PDF in the absence of smoothing obtained both 
with the formula (27) and with the reconstruction formula (19). As in ID, the reconstruction method (that 
neglects the a-dependence of the S p parameters) leads to a less extended high density (exponential) tail, but 
does not change very much the low density behavior. When the filtering effects are taken into account the 
density PDF, however, tends to have weaker non-Gaussian features both at the high- and low-density tail. 
The recipe to build density PDF in ZA with smoothing for arbitrary n is consistent in substituting eqs.(34), 
(31) into the general reconstruction formula (19), but valid for small a only. 

4. SUMMARIZING PERTURBATION SERIES FROM EXACT DYNAMICS 

In the previous section we calculated the statistics using the Zel'dovich approximation. Beyond the 
Zel'dovich approximation, that is when one assumes only the single stream regime, the general form of 
the density PDF, that would be the counterpart of (27) in ZA, has never been obtained. However, using 
perturbation theory for Gaussian initial conditions, it is possible to derive cumulants of the density PDF. 
The approximation which is admitted to apply the perturbation theory, is that the gravitational clustering 
at sufficiently large scales can be considered in the single stream regime ignoring small scale details. 
This approach to derive cumulants in quasi-linear dynamics has been intensively used in the literature 
(Peebles 1980, Fry 1984, Goroff ct al. 1986, Bouchet et al. 1992, Juszkiewicz et al. 1993, etc.). Within 
this approximation to the complicated actual dynamics, the closed form for the generating function of the 
cumulants is rigorously derived in the limit of small a (Bernardeau 1992, 1994). Since the cumulants very 
slowly increase with a while a < 1, at least for n w -1, these results can be extrapolated across the whole 
quasilinear regime, which makes them very useful. In this section we first recall the results of calculation 
of S p (0) and the PDF based on the perturbation theory without final smoothing, as it was derived in 
Bernardeau (1992). Then we report the corresponding results when final smoothing is included (Bernardeau 
1994), which are the most important for practical applications. Then we show how the general method of 
the Edgeworth expansion can be used with the cosmological dynamical perturbation theory. 

4.1. Reconstruction of PDF through S p in the small a limit from the exact dynamics 



(35) 



(36) 
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In ID, the large-scale limit (e.g. [15]) of the parameters S p is exact. In 3D the Zel'dovich solution is 
really an approximation, at any stage of the dynamics, so that the large-scale limit of the parameters <Sjf (0) 
from the ZA is also an approximation to the exact S p (0). However, it is possible to calculate the values of 
S p (0) from the exact single stream cosmological dynamics, without having to know the shape of the density 
PDF a priori. 

The principle of the calculation has been given by Bernardeau (1992) and is summarized in Appendix 
C. In this section we recall the results that have been obtained to compare them with those obtained from 
other methods. One can use the basic single stream cosmological equation of gravitational instability, and 
seek the solution in the form of perturbation series, e.g. 8 — J2"^Li see Appendix C. Then there is a 
set of basic equations in each order p for 8^ p \ Let us define the following connected (normalized) averages: 
Up = (8^(S^) p ) c /a 2p . Next, we construct the generating function of these averages 

oc 

The choice of that particular combination is based on the observation that there is a closed equation for 
Q{t). Indeed, multiplying equations for by (8^) p and averaging the result, and then summarizing the 
hierarchy of equations, one can derive the single equation for Q(t) (eq. [C7]). It is remarkable that this 
equation does not contain space derivatives, and has exactly the same form as the equation of the spherically 
symmetric collapse of the "overdense" Q(t) with a "scalar factor" r. A good approximate analytical solution 
of the "spherical collapse" is given by 

^" (l + r/W- 6 " 1 - (38) 

The function Q(t), in principle, depends on the cosmological parameters, but weakly. The closed form (38) 
is actually a fitting formula that turns out to be extremely accurate for any cosmological parameters. 

The generating function <p(y) of the S p (0) parameters (18) can then be built with the generating function 
(38) of Up parameters. In Field Theory the averages like (8^ (8^) p ) c are known as vertices, the averages 
like the cumulants (8^) c then are trees. They can be represented by diagrams where factor a corresponds 
to the lines, and factor u p to the vertices (for a vertex connecting p lines, sec Appendix C for details). There 
is a very useful and deep general result which links the generating functions of vertices and cumulants: the 
generating function of cumulants is connected to the generating function of vertices through the Legendrc 
transformation (see Bernardeau & Schaeffer 1992 for references). In the cosmological context the Legendrc 
transformation reads as 

^{y) = y + yQ{r) + \r 2 , r = -yg'(r). (39) 

Now using (39) and (20), one can derive the whole series of the S p (0) parameters. For instance, we have 

34 

S 3 (0) = -«4.9, 

60712 V°) 
S 4 ( ) = 2££f « 45.9. 
v ; 1323 



This approach allows to derive the parameters S p for a = 0. Direct perturbative calculations also allow 
to obtain the parameters £3(0), £4(0) (Peebles 1980, Fry 1984) and, additionally, should admit derivation of 
their further a dependence. Unfortunately, no a corrections of S p (a) were obtained from the perturbation 
series so far, but analysis of the numerical simulations indicates that a dependence of S p (a) is rather weak in 
the quasilinear regime (Juszkiewicz et al 1993, and hereafter §5), at least for the interesting cases n = 0, — 1. 
Outside of the quasilinear regime, S p (a) dependence might be noticeable, as some numerical simulations 
indicate (Lucchin et al. 1994). Note, however, possible impacts of the numerical effects on these calculations 
(Colombi et al. 1993). 



We can reconstruct the density PDF by substituting (39) in the general formula (19). The resulting 
shape of the density PDF is presented in Fig. 7. Note that this PDF depends on a only and not on the 
power spectrum, as it was in the ZA without smoothing. 



4.2. The final smoothing in S p and PDF for small a 

It is also possible to get the full series of the (0) parameters when the final smoothing, with a top-hat 
window function, is taken into account. We can do it via the same method we used earlier for the ZA in 
Sec. 3.6. The generating function of the cumulants is defined via a 5 s -function given by the relationship, 

where Rq is the filtering scale and a(Ro) is the rms density fluctuation at this scale. The generating function 
of the cumulants is then given by 

V s {y) = y + yG s {T) + l -T\ T = -y^-g s (r). (42) 



For instance the coefficients S§(0) and Sf (0), for a power law spectrum of index n, are 

34 

y 

60712 62, „ s 7 



Sf(0) = y-(n + 3), 



^ (0) -T32f-f ( " + 3)+ 3 (n + 3)2 - 



(43) 



We derive (43) from the generating function (41). The same coefficients S§ (0) and Sf (0) as function of 
n were derived independently directly from the perturbation theory (Juszkiewicz et al. 1993, Bernardeau 
1993), which confirms the validity of the general formula (41). 

Substituting (42) into (19), we can reconstruct the density PDF under all three assumptions. Now it 
also depends on the power spectrum. The properties of the resulting density PDF are given in details by 
Bernardeau (1994), and are recalled in the Table of Sec. 6 below. 

Despite the fact that the a dependence on the S p parameters has been neglected, it seems that the 
formula (19), with the generating function from (42), is the most reliable analytical expression for the 
density PDF in the mildly non-linear regime (see Sec. 5) for n >, — 2. 



4.3. The Edgcworth expansion in 3D case 

The accuracy of the Edgeworth decomposition in a realistic case of 3D exact single stream dynamics 
with final smoothing is worth examining. Using the expressions of 5*3 and S4 (eq. [43]) we can use the 
decomposition (22) up to the second order in u. Note that the use of this decomposition up to the third 
order would require the determination of the first er-correction ~ er 2 in S3 which we do not know yet. 

In Fig. 8 we plot the density PDF in 3D case reconstructed from the Edgeworth expansion (22) with 
S p (0) from (43) for n = -1, and compare it with the PDF from (42) and (19), for a = 0.3 and a = 0.5. We 
can see that a couple of iterations of the expansion (22) reproduce the peak of P(6) in the interval \5\ ^ 0.5 
around it for small a relatively well. It reproduces well the shift of the maximum towards the low density. 
It fails to reproduce P(S) outside of this interval. For a given value of a, each next o iteration improves the 
approximation quite slowly. The reconstruction is rapidly worsening as a increases, and in practice is useless 
for a > 0.5. 

As the S3 and S4 parameters are lower than in ID, the accuracy of the decomposition in 3D case is 
better. Correspondingly, in 3D case with the final smoothing, the accuracy of the Edgeworth decomposition 
is also better for smaller 53 and S4, i.e. for larger index n. Additionally, from the comparison of ID and 
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3D cases, one can infer that the the Edgeworth decomposition is of interest when the actual skewness of the 
distribution function is small, that is when £3 o is small. 



5. COMPARISON WITH PDF AND S p FROM N-BODY SIMULATIONS 

We have derived analytically two forms (27) and (19) of density PDF in two reasonable dynamical 
approximations, the ZA and the perturbation theory in a single stream regime. Apparently, there is no 
universal formula for the PDF in general. Now we compare our analytical results with those from cosmological 
N-body simulations. The distribution function (27) has already been tested in a previous work, Kofman et 
al. (1994) with a density field that had been filtered with a Gaussian window function. From the theoretical 
point of view, the effect of filtering is better known for a top-hat window function. We therefore run a new 
series of tests with the top-hat smoothing. 

We used a large numerical simulation kindly provided by Couchman (Couchman 1991). The simulation 
has been made in a box of 200 h~ 1 Mpc size with periodic boundary conditions, and contains 2.1 x 10 6 
particles. It used an adaptive P 3 M code and the initial conditions correspond to a CDM power spectrum 
with fl = 1, Hq = 50km s _1 Mpc -1 and the bias parameter b « 1.0. Wc made three filterings with a top-hat 
window function at two different time steps, at redshifts z = 0.6 and z = 0. The three different filtering 
radii we choose were, 5 /i _1 Mpc, 10 /i _1 Mpc, and 15 h~ 1 Mpc. The errors for all the measures that have 
been made have been determined by dividing the simulation box in eight equal subsamples and by making 
eight different measurements. 

5.1. Calculation of S p from N-body simulations 

The first test is the determination of the parameters S3 and S4 compared to the theoretical predictions. 
We calculated them from the counts of particles in the ensemble of 50 3 spheres disposed on a grid. Thus, 
it corresponds to a spherical top-hat filtering. As the number of particles is quite significant, the shot 
noise effects are negligible and have been neglected to compute the moments of the measured distributions. 
The resulting values of S3 and S4 are plotted in Fig. 9 as a function of a. For each filtering radius, 
we calculated the initial effective index n of the power spectrum to derive the expected value of of S p (0) 
coefficients from formula (43). For the three filtering radii 5 ft. _1 Mpc, 10 ft. -1 Mpc, and 15 ft. _1 Mpc we get 
correspondingly n ~ —1.3,-1.0,-0.7. Correspondingly, three values S^O) « 3.2,2.9,2.6 and three values 
64(0) ~ 17.7, 13.3, 10.6 at a = plotted in Fig. 9 without error bars are these theoretical predictions. 

Each curve is related to a different filtering radius: circles, squares and triangles correspond respectively 
to the smoothing with 5 /i Mpc, 10 /i Mpc, and 15 /i Mpc. First point without error bars on each curve 
is the theoretical prediction at er = we just discussed, which also can be interpreted as the initial time 
step at redshift z — > 00. Two other points on each curve correspond to two other time steps, at the redshift 
z = 0.6 and at present z = 0. It can be seen that the theoretical prediction (43) - how S p (0, n) depend on 
the initial index n - is well reproduced by the extrapolation of the numerical curves backward to z — > 00, 
i.e. to a — > 0. Moreover, for a given smoothing radius, the S p parameters do not exhibit any a dependence 
within the error bars. This result makes the theoretical derivation of the S p (0) series particularly attractive. 
It is also obvious that error bars are bigger for bigger filtering scale. Other numerical results by Bouchet & 
Hernquist (1992), and Lucchin et al. (1993) indicate, however, that a variation with a may be significative, 
especially in the nonlinear regime, and for low values of n. 

5.2. Calculation of PDF from N-body simulations 

The second numerical test is the construction of the density PDF from the N-body simulation. We used 
the top-hat filtering of particle distribution at different time steps z = 0.6 and z — for different filtering 
radii 5 /i _1 Mpc and 15 h~ 1 Mpc. Top-hat smoothing allows us to construct the PDF as count-in-cell statistics 
for the spherical geometry of cells. Fig. 10 shows the Eulerian density PDFs at z = 0.6 for two smoothing 
radii 5 h~ 1 Mpc and 15 h~ 1 Mpc, which corresponds to a = 0.92 and a = 0.29; as well as the PDFs at z = 
for the same smoothing but with a = 1.52 and a = 0.47. This choice of parameters covers a broad range 
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of non-linear stages, 0.3 ^ a ^ 1.5. The error bars are the standard deviation of the mean in eight equal 
subsamples. We plot P(g) over wide range of density up to g = 10. 

We compare the numerical PDF with the two most advanced analytical predictions. First, we plot the 
density PDF from (27) derived in ZA without smoothing, for four corresponding values of a. Formula (27) 
approximates numerical PDF with top-hat smoothing very well up to a ^ 0.5 in the range < g 3, then 
starts to overestimate the high-density tail. This is the range of validity of the underlying assumptions of the 
Zel'dovich approximation without final smoothing. For larger a the approximation (27) is slightly worsening, 
and is out of applicability in the multiple stream regime a > 1. Our conclusion confirms of that of Kofman 
et al. (1994), based on the Gaussian filtering. 

We also plot the density PDF from the exact perturbation theory in the single stream regime and with 
the final smoothing. Assuming that S p parameters are constant in quasilinear regime, we substitute the 
generating function f S (y) from equation (42) into the reconstruction formula (19), and calculate P(g) for 
corresponding values of a and n. The results are presented in Fig. 10, and show a remarkable agreement 
with the numerical density PDF over the whole range of g. We extrapolate the theoretical PDF from the 
perturbation theory for non-small a beyond the range of its validity. However, the agreement with the 
numerical PDF is striking up to the maximal used values a ~ 1.5 and g w 10, and so far there are no signs 
of deviation even for higher a\ Plausible explanation of why formulae (19), (42) work so well is that S p (a, n) 
parameters depends very weakly on a up to moderate er, at least for n ^ — 2. It has been explicitly checked 
for S$ and S4 but it should also be true for any of them otherwise the low- and high-density tails would 
not have been reproduced so accuratly. Then it is not too surprising that the density PDF (27) from ZA 
for which the S p parameters are not constant (see Fig. 6) does not fit well the low- and high-density tails 
for n £ — 2. For n ^ — 2, however, the S p (a) parameters may have a stronger dependence on er, and the 
ZA could then provide a more reliable PDF. In any case it would be interesting to check these theoretical 
predictions against the PDF from N-body simulations for higher a. 

5.3. Fitting by the log-normal distribution 

As it was noted a long time ago by Hubble (1934), the galaxy count distribution in the plane cells 
on the sky might be well described by the log-normal distribution. The log-normal distribution fits the 
observed galaxy PDF from 3D surveys as well (Hamilton 1985, Bouchet et al. 1993, Kofman et al. 1994). 
The log-normal density distribution reads as 

^, o? = ln(l + " 2 )- (44) 
Q 

Kofman et al. (1994) found that the log-normal distribution is an excellent approximation to the density 
PDF from N-body CDM simulation for moderate values of er in the used range q < 5. In Fig. 10 we also 
compared the density PDF from N-body CDM simulation with the log-normal distribution. We also found 
a striking agreement between the log-normal PDF and that from N-body CDM simulation for the tested 
values of 0.3 < er < 1.5 in the tested range of g < 10! 

Such a remarkable fitting inspires the thought that there might be a strong dynamical reason to manifest 
the log-normal features of the density PDF. For instance, Coles & Jones (1991) argued for the log-normal 
mapping of the linear density field to describe its non-linear evolution. Their log-normal model is universal 
for any spectral index n. Unfortunately the log- normal mapping does not work (Coles et al. 1993). 

Why does the log-normal density PDF work so well? 

We argue that the log-normal successful fit can be seen as a mere coincidence due to the shape of the 
CDM power spectrum at moderate a. The log-normal PDF is not a universal form of the cosmic density 
PDF due to the non-linear dynamics, but is rather a convenient fit for the particular region in the plane of 
(er, n)-paramctcrs. This region of (er, n)-paramctcrs includes the CDM model at moderate er. It explains why 
the log-normal PDF fits the results of N-body CDM simulations. Consequently, the "log-normalish" features 
of the observed density PDf mean that the realistic cosmological model corresponds to that (er, n)-region, 
i.e. close to the CDM model in this respect. 
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Let us consider the properties of the cumulants of the log-normal distribution. The first two S p 
parameters of this distribution as a function of a are 

^ og (a) =3 + a 2 , (45a) 

S%*(a) = 16 + 15a 2 + 6a 4 + a 6 , (456) 

for an arbitrary a. For the sake of completness note that the log-normal distribution has S p ° e (0) = p p ~ 2 , 
and its t/ log (r)-function is simply cxp(— r). 

Parameters S l ^ s (a) and S l 4 S (a) are plotted on Fig. 11 as dashed lines. We compare curves (45) with 
the values of S 3 and 5 4 from N-body CDM simulations for four values of a = 0.29,0.47,0.92,1.52. The 
log-normal curves (45) are shown to match the CDM parameters for moderate values of a reasonably well. 
The Edgeworth expansion (21) then helps to explain why the shape of the two density PDFs - from N-body 
CDM simulation and from (44) - are in good agreement if there is good agreement between S p parameters. 
Consequently, we predict the agreement worsening as far as S p -s are departing from each other. The accuracy 
of the log-normal distribution is good essentially due to the particular slope of the CDM power spectrum at 
corresponding smoothing scales at moderate a. A less steep power spectrum, for instance, would not have 
led to the same level of agreement. The formula (42) with (19) shows what the expected dependence of the 
density PDF with n is (see Bernardeau 1994 for more details). The log-normal distribution is expected to fail 
for the higher values of a where the density PDFs from N-body simulations exhibit a power law behavior (see 
for instance, Bouchet & Hernquist 1993), which is not the case for the form (44). The low- and high-density 
tails of the log-normal distribution are very different from those given by the reconstruction formula (19), 
as it is shown in the Table below. 

In practice, however, for an observed power spectrum n « — 1 at moderate a, the log- normal distribution 
is a very effective and simple fit for the density PDF in the mildly non-linear regime. In accompanying paper 
(Kofman & Bernardeau, 1994) we address the domain of validity of the log-normal fit more specifically. 

6. CONCLUSION 

6.1. Theoretical framework 

Let us recall the assumptions that have been made throughout the paper to derive the various presented 
results. As we noted in Section 3.1, the actual gravitational dynamics for the realistic cosmological scenarios 
is quite complicated and includes the superposition of hierarchical pancaking and clustering across the vast 
range of cosmological scales. In the non-linear regime, the basic equations in terms of the continuous cosmic 
fileds are multiple stream. However, at large scales, ignoring small scale substructures, we expect that the 
gravitational clustering is simpler, and might be approximated by the single stream regime. Unfortunately, 
this transition to the large-scale single stream description (i.e. the demonstration that, indeed, the small 
scale details can be ignored for the large-scale dynamics) has never been done, even in the linear regime. The 
common (but formally unjustified yet) belief is that after the large-scale filtering the perturbation theory, 
or truncated ZA can be applied. It relates to a lesser extent to the N-body simulations with inevitable 
truncation of the genuine power spectrum. We are also working within this assumption, basically because 
the results we derive are in good agreement with the N-body simulations. However, we clearly understand 
the nature of the approximation and the need to justify it. 

The usual practice is to filter the initial fluctuations in order to model the large scale dynamics. Thus, 
our first assumption is that 

la) the initial fluctuations are smoothed to ensure being in the single stream regime, the particular 
filtering scale R is controlled by the rms density contrast, a = a(R). 

We used two theoretical approaches (in particular, to derive the density PDF). One is to assume that 

2a) the truncated Zel'dovich approximation can be used to describe the large-scale dynamics. 

The results which might be derived under these two assumptions are given in §3. Note that the ZA can 
be used whatever the nature of the initial conditions is, Gaussian or not. 
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The other approach is the perturbation theory for the cosmological equation of continuity, the Euler and 
the Poisson equations. This method intrinsically relies on the hypothesis of Gaussian initial conditions. But 
in any case, even in the single stream approximation, the equations arc fully non-linear and it is difficult, if 
not impossible, to find an exact solution. For instance, the general form of the S p (a) parameters as functions 
of g are beyond our current skills of analytic calculations. We are then led to derive them in their leading 
order in cr, that is 

2b) the cumulants are derived in the limit of small cr. 

This calculation can be done rigorously from the summarized perturbation series. Note, that in the 
paper Bernardeau (1992) the PDF was derived under the assumption la) and 2b), contrary to how this 
paper is sometimes quoted. 

So, within the assumption la) and 2a) or la) and 2b) it is possible to make interesting theoretical 
predictions, for instance, get the dynamical derivations for the one-point density PDF and moments. 
However, these calculations do not take into account the final filtering that, in practice, cannot be avoided. 
This is a crucial step since it alters the statistical properties of the cosmic fields. If the scale of the final 
filtering is larger than the scale of the initial filtering, then the initial filtering is irrelevant, and we can 
replace the assumption la) by the assumption that 

lb) this is the final smoothing only, which assures that the large-scale dynamics can be accurately 
described by the single stream approximation. 

Even within the simplified dynamics given by the ZA it is no longer possible to derive the density 
PDF under this last assumption. The only known approach that allows to do this comes from perturbative 
calculations, since it makes it possible to derive the leading order of the cumulants of the final density field. 
It is then a considerable improvement compared to previous results. The assumptions lb) and 2b) then lead 
to a density PDF in very good agreement with the numerical results. 

6.2. Results 

We derived the density PDF, cumulants in forms of S p (cr)-parameters and their generating functions in 
case of ID gravitational dynamics, in the 3D Zel'dovich approximation with and without final smoothing, in 
the perturbation theory extrapolated over mildly non-linear regime, with and without final smoothing, for 
the log-normal distribution, and from cosmological N-body CDM simulation. We summarize the quantitative 
results in the Table. 

We also present new qualitative results stemming from our study. As one can see from the Table, the 
values S p (0) are constants which characterize the particular dynamical model. For example, these values 
are different for ID gravitational instability, 3D Zel'dovich approximation, single stream 3D gravitational 
instability, sea waves dynamics, etc. The common trend is that S p is rapidly increasing with p, but at a 
different rate. The cr-dependence of the S p (a) parameters also characterizes the particular dynamics, and 
as we illustrated for different models, might also be quite different. The smoothing effects introduces an 
extra dependence on the power spectrum index n, S p (<r,n). We can expect S p (<r, n) from different models 
can coincide in some regions of parameters (cr, n) . It allows to construct simple fitting formula, as we found 
it for the log-normal distribution. 

Sp(0)-coefficients are related to the particular non-Gaussian statistics which emerges from the non- linear 
dynamics. The usual practice of their direct measurements from observations is significantly affected by the 
final sample volume. We, however, learned from the Edgeworth expansion that S p (0) coefficients are related 
to the shape of the PDF peak. It gives us an alternative method to evaluate the skewness and kurtosis 
by measuring the PDF around its maximum, which is statistically more robust. This approach might be 
interesting in other contexts, such as to constrain the skewness of the cosmological AT/T fluctuations, or 
skewness and kurtosis of the cosmological velocity field and its divergence. 

There is a deep connection between the generating function of S p (0) parameters and ^(r)-function 
- the generating function of vertices, and the non-linear dynamics of the spherical overdense fluctuation 
in the particular dynamical model. The solution of the last problem in the particular dynamical model 
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gives the ^(r)-function, which allows to reconstruct the density PDF in the limit of small a. This is a 
remarkably simple and general prescription. From the Table note a simple generalization of the G(t) for 
the models without final smoothing: Q{t) = (l + r/a) - " — 1. For the Zel'dovich approximation in the 
space of TV-dimensions, a = N . For the 3D perturbation theory a ~ 1.5. In the limit W > 1 we get 
Q{t) = cxp(— t) — 1, which coincides with that of the log-normal distribution. 

We have derived the cosmological density PDF in two different approximations. The one based on ZA 
gives the density PDF in the early non-linear regime, a £ -5, and is expected to improve for the models with 
n ^ — 2, for which pancaking is more pronounced. This approach remains irreplaceable for models with 
non-Gaussian initial density fluctuations. However, we found that the best theoretical model for the density 
PDF evolving from Gaussian distribution, in the case of n }t — 2, is based on the perturbation theory when 
the final smoothing is included. This PDF works remarkably well for significant range of a in the mildly 
non-linear regime. Both approaches provide us with an efficient machinery to deal with one-point statistics 
for a broad range of models. In the Table we summarize the properties of the various approximations that 
have been presented and used in this paper. It shows, in particular, that the low-density tail is not affected 
by the reconstruction method. On the other hand, the shape of the high-density tail is dramatically modified. 
Both tails are also slightly modified, where the smoothing effects are taken into account. 

It was noted earlier and further confirmed that the log-normal distribution is an excellent fit for the 
density PDF from CDM non-linear dynamics. We found an explanation of this mystery, based on the 
properties of cumulants. The log- normal distribution fits well in the particular range of the parameter space 
(n, a) around n ~ —1, a ~ 1/2, and worsening outside of this region. By chance, the popular CDM model 
at moderate a corresponds to this region. Thus, some "log-normal" features of the observed density PDF 
would mean that realistic cosmological model is close to that range of parameters. 
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APPENDIX A: Statistical technics 

In this appendix, we give technical definitions of interest for statistical studies, such as the moments 
of the density distribution function, its cumulants, the generating functions of the moments and of the 
cumulants. Relationship between those quantities are also given, -lz The moments (8 P ) of the distribution 
function, P(p), are given by the integrals, 

/>oo 

(F) = / dgP(g)(g-ir. (Al) 
Jo 

The cumulants can then be obtained from the moments. The p th cumulants, (S p ) c , is defined recursively 
from the p th moment following the rules, 
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In general, to obtain the p th cumulant, one should consider all the decompositions of a set of p points in its 
subsets (but the one being only the set itself); multiply, for each decomposition, the cumulants corresponding 
to each subset and subtract the results of all these products obtained in that way from the p th moment. 

Each cumulant, in a given order, actually contains a piece of information about the shape of the PDF 
that cannot be derived from the lower order cumulants. For instance, the second cumulant gives the width 
of the distribution, the third cumulant measures its asymmetry, the fourth - its flatness, etc. The Gaussian 
distribution is characterized by only one non-zero cumulant (whereas all even moments are nonzero): the 
second one, which gives its variance. However, in general, a distribution function will be characterized by the 
whole series of cumulants. Therefore, to treate more easily the series of the cumulants, we arc led to define 
other mathematical tools of great practical interest. There are the generating functions of the moments and 
the cumulants. The generating function of the moments is defined by, 

oo p 

The generating function of the cumulants, C(A), is defined in a similar way, 



p=2 y 

One fundamental result of the statistics is that M.{jjl) and C(p) are connected in a simple way. Indeed, 
we have 

MQj) = cx P [C(m)]. {Ah) 

We omit here the rigorous demonstration of this property (see, e.g. Balian & Schaeffer 1989). By expanding 
cxp[C(/x)] with respect to fi, one can easily verify that the first few moments are given correctly. 

Using this property, it is then possible to relate the generating function of the moments, or of the 
cumulants, to the shape of the PDF. Indeed, we have 
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APPENDIX B: Derivation of the Edgeworth expansion 

In this appendix we present the derivation of the Edgeworth expansion. It is based on the reconstruction 
(19) for the density PDF. The generating function ip(y) can be expanded with respect to y (eq. [20]). We 
then have to expand the non-Gaussian part of the exponent in (19), 
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This expansion should be made with respect to both y and a, assuming they are of the same order. The 
resulting value of the density PDF requires the determination of integrals of the form 
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where H n {v) are the Hermite polynoms 
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The resulting form of the density PDF is 
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APPENDIX C: S p (0) from the Perturbation Theory 

In this appendix we give a brief sketch of the derivation of some technical details used Sec. 4.1. 
Cosmological gravitational instability of the perfect fluid without pressure in the single stream regime is 
described by the continuity equation 



the Euler equation 



and the Poisson equation 



| + 3^ + iv x (pv)=0, 

dv 1 1 

^7 + -(vVx)v + Hv = — V x </>, 

at a a 
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(Clb) 
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It is convenient to use a(t) as a new time variable, and the velocity potential $ for the potential flow under 
consideration, see §3.1 for definition. Let us consider the case of the Einstein-de Sitter Universe. Then three 
equations (CI) can be rewritten in form of two equations: 

d 

a— <5(x, a) + (1 + 5(x, a))A$(x, a) + V<5(x, a).V$ = (C2a) 

and 

d 1 
a— A$(x, a) + -A$(x, a) + V$(x, a).V(A$(x, a)) 

3 3 (C2b) 

+ 2^ *,a/3(x, a)$ Q/3 (x, a) + -<5(x, a) = 0, 

a, 0=1 

with $ jQ/ g(x, a) = (1/aH) 2 9 2 $(x, a)/9x ct 9x / 3 where x Q is the component a of x. 

Let us seek the solution of the equations (C2) in the form of perturbation series, e.g. 8 — Y^Li 
V 2 $ = X)^Li(^ 72< ^)' P ^: etc - Then we obtain the hierarchy of equations for values S (p \ (V 2 $)( p ) etc., in each 
order p of the perturbation series. Then, at the lowest order of a, the cumulant ($ p ) c is given by a term in 
the form 

<*•>„= E (h* (p[l]) ) c (G3a) 

combinations, p(i) i—1 

where the sum is taken over all the possible combinations p(i), i = 1 . . .p, for which 

p 

i=l 

This particular result is correct due to the hypothesis of the Gaussian initial conditions. The terms in (C3b) 
are all products of the rms density fluctuation at the power 2p — 2 by some combination of the vertices, 

Vp = (5^\5^) n )Ja 2p . {CA) 

(and respectively \x v for V 2 <3>, etc.) The parameters S p (eq. [7]) for a — are then given by a combination 
of the vertices, v p , that technically corresponds to a tree summation. The relationship between the S p 
parameters and the vertices can be written in a closed form at the level of their corresponding generating 
functions. Let us define the function Q{t) by 
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(a similar one can be denned for other vertices). Then we have 

l p(y,0)=y + y g(r) + ^-, r = -y Ag( T ), (CQa) 

where <p(y,a) is defined in (20). 

The problem then is to derive the function Q(t) from the equations (C2). F.B. 1992 derived useful 
properties of the generating functions of vertices. Multiplying equations (C2a) by (8^) p and averaging the 
result, and repeating similar operation with eq. (C2b), and then using properties of the generating functions 
of vertices, one can finally get the single equation for G(t): 

d 2 „ 4 / d „\ 2 3,. _ d „ 3 



with Q{t) ~ — t when r — > 0. 

It is remarkable, that this equation does not contain any space derivatives, and has exactly the same form 
as the equation of the spherically symmetric collapse of the "overdense" Q{t) with a "scalar factor" r. The 
analytical solution of the equation (C7) is well known. When r < 

HISS- —!Gc-'>f 



and when r > 



^ 9(sinh<9-60 2 , 3 /3, . , „ „.\ 2/3 



It turns out that the form (38) for <7(r) w 1/(1 + r/1.5) 15 — 1 is a very good fit to the exact solution 
(C8). 

In general case of arbitrary background cosmology the function Q(t) depends on the cosmological 
parameters. Fortunately, this dependence is weak, and the form (38) can be accurately used in general 
case. 
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Table and Figure captions 

Table 1. Properties of the density PDF obtained with the various approximations described in the 
text. The first two columns give the values of the S3 and 5*4 parameters. In general case they depend on a. 
A couple of first terms of the a-expansion are given in cases they are known. The third column gives the 
shape of the function Q(t) which is used to derive the function f(y), involved in the reconstruction formula 
(19). The last two columns give the high- and low-density asymptotas of the resulting distributions. 

Fig 1. Shape of the density PDF for ID dynamics and for a = 0.5 (eq. [4]) (solid line). The dotted, 
dashed, long dashed and dotted dashed lines show the shape of the density PDF defined in Eq. (6) for 
respectively g c =5, 7.5, 10, 12.5. 

Fig 2. coefficients after regularization for the ID dynamics. The filled circles correspond to the 
theoretical limits for a — ► (Eqs. [15]), and the thick solid lines correspond to the theoretical curves (Eqs. 
[16]). The dotted, dashed, long dashed and dotted dashed lines are for a sharp cutoff (Eq. [5], left panel) or 
an exponential cutoff (Eq. [6], right panel) for respectively g c =5, 7.5, 10, 12.5. 

Fig 3. The shape of the density PDF as obtained from eq. (4), thin solid lines, and by the reconstruction 
method [eq.(19)], thick solid lines, for a = 0.3 (left panel) and a — 0.5 (right panel). 

Fig 4. The Edgeworth expansion (Eq. [22]) compared to the form (4) for a = 0.1 (left panel) and 
a = 0.3 (right panel). The dashed line corresponds to the case when the skewness only is taken into account 
(up to a correction), the long dashed lines when both the skewness and the kurtosis are taken into account, 
up to a 1 corrections, and the dotted line when the expansion is made up to the a 3 order (Eq. [22]). 

Fig 5. Shape of the density PDF for a = 0.3 (left panel) and a = 0.5 (right panel) for various methods 
based on the Zel'dovich approximation. The solid line is the shape obtained by the Zel'dovich approximation, 
(Eq. [27]); the other curves have been obtained by the reconstruction formula (19) by assuming that the 

ratios (d p ) c /{3 2 ) P = Sp equal their low-er limit in the ZA without taking into account the smoothing 
effects (dashed lines), and taking into account the smoothing effects for n = — 1 (long dashed lines). 

Fig 6. parameters after regularization for the 3D dynamics. The dotted, dashed, long dashed and dotted 
dashed lines are for the sharp cutoff (left panel) or the exponential cutoff (right panel) for respectively p c =5, 
7.5, 10, 12.5 in the regularized expression (27). The filled circles correspond to their theoretical limit at small 
a and the squares to their theoretical limits when the filtering effects are taken into account, for n = —1. 

Fig 7. Shape of the density PDF for a = 0.3 (left panel) and a = 0.5 (right panel) from the 
reconstruction formula (19) using the value of the cumulants from (38) for the exact dynamics in the a — > 
limit. The solid line is the density PDF obtained by assuming that the values S p equal their low a limit 
without taking into account the smoothing effects. The dashed line is obtained when the smoothing effects 
are taken into account assuming that the initial power spectrum is P(k) oc fc _1 . The long dashed line is the 
log-normal distribution. 

Fig 8. The Edgeworth expansion (Eq. [21]) compared to the density PDF obtained with (42) and 
the reconstruction formula (19) for n = —1, a = 0.3 (left panel) and a — 0.5 (right panel). The dashed 
line corresponds to the case where the skewness only is taken into account (up to a correction), and the 
long dashed lines to the case where both the skewness and the kurtosis are taken into account, up to a 2 
corrections. The long dashed line is the lognormal distribution. 
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Fig 9. The coefficients 83(a) and 84(a) as functions of a in a CDM numerical simulation. The circles, 
squares and triangles correspond respectively to the smoothing radii of 15, 10, 5 ft _1 Mpc. Three points at 
each curve corresponds to three different time steps, at z — 00, 0.6, 0. The values at a = (or z — > 00) are 
the theoretical predictions from (43) taking into account the filtering effects. 

Fig 10. The density PDF for CDM initial conditions. The points are measured in a numerical simulation 
at two different time steps corresponding to a/a — 0.6 (upper panel) and a/a = 1 (lower panel) and at 
two different smoothing radii Rq = 5/i _1 Mpc and Rq — 15/i _1 Mpc. The rms density fluctuation are then 
respectively a = 0.92 and a = 0.29 in the upper panel and a = 1.52 and a = 0.47 in the lower panel. The 
error bars have been obtained by dividing the sample into eight subsamples. The solid line is the prediction 
given by (19) when the smoothing effects are taken into account (with Eq. [42]). The dashed line is the 
prediction (27) from the ZA and the long dashed line is the lognormal distribution (44). 

Fig 11. The measured values of the 83(a) and 84(a) parameters in the CDM simulation as functions 
of a. The solid lines correspond to the values of these coefficients from formula (43) for n = —1. The 
thick dashed lines correspond to the log- normal distribution, cq. (45). The S3 and 54 values in CDM and 
log-normal models overlap for moderate a ~ 0.5 only. 



